The 3D Grazing Collision of Two Black Holes 
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We present results for two colliding black holes (BHs) , with angular momentum, spin, and unequal 
mass. For the first time gravitational waveforms are computed for a grazing collision from a full 
3D numerical evolution. The collision can be followed through the merger to form a single BH, and 
through part of the ringdown period of the final BH. The apparent horizon is tracked and studied, 
and physical parameters, such as the mass of the final BH, are computed. The total energy radiated 
in gravitational waves is shown to be consistent with the total mass of the spacetime and the final 
BH mass. The implication of these simulations for gravitational wave astronomy is discussed. 
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The collision of two black holes (BHs) is considered by 
many researchers to be a primary candidate for generat- 
ing detectable gravitational waves, and hence is the focus 
of attention for many research groups worldwide. As the 
first generation of gravitational wave detectors J|] , with 
enough sensitivity to potentially detect waves, is com- 
ing online for the first time within a year, the urgency 
of providing theoretical information needed not only to 
interpret, but also to detect the waves, is very high. How- 
ever, even in axisymmetry, the problem has proven to be 
extremely difficult, requiring nearly 20 years to solve in 
even limited cases (e.g. Q-^]). In full 3D progress has 
been rather slow due to many factors, including (but 
not limited to) unexpected numerical instabilities, lim- 
ited computer power, and the difficulties of dealing with 
spacetime singularities inside BHs. The first true 3D sim- 
ulation of spinning and moving BHs was performed in . 
In Jt|, the two BHs start out very close to each other, 
closer than the separation for the last stable orbit, and 
the evolution proceeds through parts of the plunge and 
ring-down phase of a "grazing collision" within a very 
short time interval. The spacetime singularities are dealt 
with by a particular choice of coordinates, singularity 
avoiding slicing and vanishing shift. 

BH excision (8|^] has allowed improvements in the 
treatment of the spacetime singularities to the extent 
that highly accurate simulations of single BHs can be 
carried out 10-14 and recent applications to the grazing 
collision of BHs show promise jl5| . One of the key limit- 
ing factors in the existing two approaches to the grazing 
collision is the achievable evolution time for which useful 
numerical data can be obtained, which due to numerical 
problems has been limited to 7M in JtJ, and to about 
9M-15M in p5[ . Here time is measured in units of the 
total "ADM" mass M of the system as opposed to using 
the bare mass m of one of the BHs. Note that the grazing 
collision of black holes without spin has also been studied 



in the close limit approximation [|16|- 

In this paper we consider singularity avoiding slicing. 
We combine the application of a series of recently de- 
veloped and tested physics analysis tools and techniques 
with significant progress made in overcoming the prob- 
lems mentioned above, including the use of much more 
stable formulations and much greater computer power. 
Early, preliminary results from this series of simulations 
have been presented in |j^,^0), but we now provide the 
first detailed physics analysis not possible previously. 

We compute BH initial data of the puncture type , 
corresponding to two BHs in orbit about each other, with 
unequal masses, linear momentum, and individual spins 
on each BH. The construction of such data sets, which 
involves solving the non-linear elliptic Hamiltonian con- 
straint equation numerically, is described in ]2l] |. A de- 
tailed survey of a sequence of such data sets including 
various physical properties is discussed in p2| . In this 
paper we choose punctures for each BH on the y-axis at 
±1.5m, masses mi = 1.5m and ni2 = m, linear momenta 
P h2 = (±2,0,0)m, and spins Si = (-1/2, 0, -l/2)m 2 
and S*2 = (0, 1, — l)m 2 . Note that the linear momentum 
is perpendicular to the line connecting the BHs, equal 
but opposite for a vanishing net linear momentum, and 
that the spins are somewhat arbitrarily chosen to obtain 
a general configuration. 

For this case, an asymptotic estimate for the initial 
ADM mass is M — 3.22m. Solving the Hamiltonian con- 
straint leads to a larger value than the Brill-Lindquist 
mass of mi + rri2 = 2.5m. The angular momentum for 
puncture data is given by (independent of the solution to 
the Hamiltonian constraint) J = 2d\ xPi+Si+S?, where 
di is the vector from the origin to the first puncture. 
The total angular momentum is therefore J = 7.58m 2 , 
which corresponds to an angular momentum parameter 
of a/M = J/M 2 = 0.73. In this configuration the indi- 
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FIG. 1. Root-mean-square value of the Hamiltonian con- 
straint on a centered cube with outer boundary at 38m and a 
gridspacing of 0.30m, 0.24m, 0.20m. The curves are rescaled 
so that they coincide for second order convergence. 



vidual spins increase the total angular momentum, so we 
call it the "high- J" case. The following discussion refers 
exclusively to this one data point in parameter space, ex- 
cept that when discussing waveforms below we compare 
the high- J case with data where the individual spins van- 
ish (mcdium-J, M = 3.00m, J = 6.00m 2 , a/M = 0.67) 
or where S\ — » —Si and S2 —> —£2 (low-J, M — 3.07m, 
J = 4.64m 2 , a/M = 0.49). 

This initial data is evolved with evolution equations 
of the "BSSN" family [^3 24|, using the implementation 
that we developed and tested for the collapse of strong 
gravitational waves to BHs in p5p . We discuss some 
reasons why certain variable choices and certain com- 
binations of the evolution equations with the constraints 
can lead to more stable evolutions than the traditional 
"ADM" system in p6|, and we do observe a significant 
improvement in numerical stability in practice. We use 
radiative boundary conditions for the outer boundary. 
The coordinate singularities at the BH punctures are 
handled as in j?],[nj by a time independent conformal 
factor. We solve the maximal slicing condition on the 
initial slice and then use the so-called 1+log slicing for 
the lapse and vanishing shift during the evolutions. 

The computer simulations were carried out on a 3D 
cartesian grid. On a 256 processor SGI/Cray Origin 
2000 machine at NCSA we were able to run simulations 
of 387 3 , which take roughly 100GB of memory (to our 
knowledge this makes them the largest production nu- 
merical relativity simulations to date). A good balance 
between resolution in the inner region and distance to the 
outer grid boundary was achieved for a grid spacing of 
0.2m, which puts the outer boundary for a centered cube 
at about 38m or about 12M. All said, the combination 
of resolution, outer boundary location and treatment, co- 
ordinate choice, evolution system and puncture method 
for the BHs allows evolution times past 30M. The low- 
est quasi-normal mode of the ring-down phase of the fi- 
nal rapid Kerr BH has a period of about 13M (17M for 
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FIG. 2. Waveform at resolutions 0.30m, 0.24m, 0.20m. 
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FIG. 3. Mode I = m = 2 of the even Zerilli function ex- 
tracted for different radii as a function of time. A wave that 
develops after the BHs collide is propagating out. 



Schwarzschild), therefore evolution times of 30M or more 
are a prerequisite for wave extraction, which was not pos- 
sible in 0,^5). The simulations do not crash at that time, 
but as we will discuss now, the numerical data becomes 
degraded due to effects of the outer boundary and due to 
grid-stretching (i.e. large metric gradients) in the vicinity 
of the BHs. Fig. [I] shows the root-mean-square value of 
the Hamiltonian constraint over the entire grid for differ- 
ent resolutions but same outer boundary location. The 
inset shows clean global second order convergence up to 
about 6M. A local analysis shows that there are large 
contributions to this average from inside the horizon, and 
that smaller errors intrude from the outer boundary, but 
the code is convergent beyond 30M. 

Since a main result of the simulations are waveforms, 
the most relevant measure and often most stringent cri- 
terion for numerical quality is convergence in the wave- 
forms. We use the gauge invariant waveform extraction 
technique, developed originally by Abrahams |^7j and ap- 
plied to the 3D case in |2S| , to extract gravitational wave 
modes of arbitrary £, m. As shown in p8| , p9[ , this tech- 
nique can be used on numerically evolved 3D distorted 
BH spacetimes to produce very accurate waveforms away 
from the BH, even if errors are rather large near the 
horizon. Here, we extract for example the nonaxisym- 
metric I — m = 2 mode, expected to be one of the most 
important modes in binary BH coalescence |3(]]. Fig. ^ 
shows for three resolutions the Zerilli function ip22 en (t) 
extracted at R — 7.8M. Up to a time t w 30M the 
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FIG. 4. A fit to the quasi normal mode determined by M 
and a shows good agreement in the frequency and decay rate 
at late times for a resolution of 0.2m. 
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FIG. 5. Even and odd wave parts showing differences de- 
pending on high, medium, and low-J data. 



dependence on resolution is rather small, which suggests 
that the resolution reaches the convergent regime. 

In Fig. [|, we show a sequence of extracted waves at 
different radii, obtained by integration over the corre- 
sponding coordinate spheres, as a function of time. The 
outermost detectors show late time problems due to spu- 
rious signals propagating in from the outer boundary, 
while the inner detectors are affected by the closeness to 
the strong field region. Note that these methods assume 
a Schwarzschild background, but they can be applied on a 
rotating BH, the primary effect being an offset depending 
on the rotation parameter a J3l[]. In Fig. |], we show the 
I = m = 2 even parity wave for the detector at R = 7.8M 
and a match to the corresponding lowest quasi-normal 
mode plus the first overtone. The values for M and a de- 
termine the quasi-normal frequency, while the amplitude 
and the offset in time are fitted. The observed period is 
13M, which is consistent with a final distorted Kerr BH 
with a/M = 0.73. Gravitational waves carry away energy 
and momentum from the BHs. For the energy, we have 

ds/dt = i/(32tt) yt=2 EL=-,(wr/*) 2 + WO 2 )- 

Integrating the / = 2, 3, 4 modes up to t = 35M, we find 
AE = 0.0323m w 1%M. 




FIG. 6. The merger of the AH. Shown are marginally 
trapped surfaces at times 2.5M, 3.7M, 5.0M, and 6.2M. The 
apparent horizon is the outermost of these surfaces. 



One of the potential insights from the detection of 
gravitational waves is the determination of the orienta- 
tion of spins in relation to the orbital motion. Fig. || 
shows the wave signature for the high, medium, and low- 
J data. For vanishing spins (medium- J), we note that 
^/>20 i s zero within numerical accuracy, while i/ , 22 d shows 
an oscillation with amplitude 3 x 10~ 5 . A comparison 
with the methods of would be useful. 

While it will be the waves that we can observe on earth 
directly, it is also interesting to compute the apparent 
horizon in the grazing BH collision. During the evolu- 
tion we use a 3D apparent horizon (AH) finder described 
in 1 32 1 to track the location of the horizon. In principle, 
the event horizon can also be located by techniques de- 
veloped in Ref. [Q, but we do not yet know whether a 
single event horizon is present on the initial slice in this 
data set. Fig. ^ shows the AH during a grazing collision. 

We compute the BH mass Mah and compare with 
the ADM mass of the initial data and the radi- 
ated energy to assess the overall energy account- 
ing. In Fig. |7[ we show the result of the calcula- 
tion of the so-called irreducible mass as a function 
of time, defined as Mj r = y/ Area ah /IQtt. The hori- 
zon mass Mah can be determined through the formula 
M\ H = (Mir) 2 + J 2 /(2M lr ) 2 , where we use J = 7.58m 2 
of the initial data. The observed upward drift in Mi r 
may be curable by excision or better coordinate condi- 
tions, but even in the present case we can estimate the 
final mass of the BH to be Mi r w 3.0m and Mah ~ 3.3m 
in this simulation. Comparing this to the initial ADM 
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FIG. 7. We show the evolution of Mi T for the high-J con- 
figuration. As the BHs merge, the area grows, and then begins 
to level, as the final BH goes into the ringdown. But numer- 
ical error associated with the grid stretching effects causes a 
spurious growth in the area, familiar in previous 2D and 3D 
studies. However, one can estimate the mass of the final BH, 
as shown by the dashed line. 



mass of the spacetime, M = 3.22m, we find consistency 
in the overall energy accounting from independent phys- 
ical measurements. The fraction of the total energy, 1%, 
that is carried away in the gravitational waves falls within 
the error estimate of this energy balance. 

In conclusion, these results indicate that for the first 
time we are indeed able to simulate the late merger stages 
of two BHs colliding, with rather general spin, mass, and 
momenta, and that we can begin to examine the fine de- 
tails of the physics. Studies of apparent horizons, wave- 
forms, and asymptotic properties show consistency in the 
analysis across strong field, near zone, and far field re- 
gions. Without more advanced techniques, such as BH 
excision, these simulations will be limited to the final 
merger phase of BH coalescence. But while that is under 
development, we can take advantage of our capabilities 
and explore this phase of the inspiral now. Our goal is 
several fold: (a) to explore new BH physics of the "final 
plunge" phase of the binary BH merger, (b) to try to 
determine some useful information relevant for gravita- 
tional wave astronomy, and (c) to provide a strong foun- 
dation of knowledge for this process that will be useful 
when more advanced techniques are fully developed. A 
new development is the Lazarus project, which provides 
an interface between full numerical relativity simulations 
and perturbative evolutions |34| , and an interface to the 
post-Newtonian inspiral phase is planned. At this time 
Lazarus allows us to move a finite time interval of full 
numerical evolution into the early stages of the merger, 
and to perform the ring-down calculation efficiently once 
the final BH enters the perturbative regime. When these 
techniques and excision are used to extend the ability of 
the community to handle the collision of two BHs starting 
from the late orbital phase, it will be important to have 
an understanding of details of the most violent merger 
phase in advance, both as a testbed to ensure that re- 
sults are correct, and because the understanding we gain 
may be useful in devising the appropriate techniques for 



longer term evolution. 
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